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We report on numerical results from an independent formalism to describe the quasi-equilibrium 
structure of nonsynchronous binary neutron stars in general relativity. This is an important indepen- 
dent test of controversial numerical hydrodynamic simulations which suggested that nonsynchronous 
neutron stars in a close binary can experience compression and even collapse prior to the last stable 
circular orbit. We show that the interior density indeed increases as irrotational binary neutron stars 
approach their last orbits for particular values of the compaction ratio. The observed compression 
, is however at a significantly reduced level. 
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I. INTRODUCTION 



The physical processes occurring during the last orbits of a neutron-star binary are a subject of intense current 
debate [|l-ll|. In part, this recent surge in interest stems from relativistic numerical hydrodynamic simulations in 



. which it has been noted [|l]-|3| that as the stars approach each other their interior density increases. Indeed, for an 
appropriate mass and equation of state, previous numerical simulations indicated that binary neutron stars might 
— , , collapse individually toward black holes many seconds prior to merger. This compression effect would have a significant 
^ ' impact on the anticipated gravity-wave signal from merging neutron stars and may also provide an energy source for 
l/") cosmological gamma-ray bursts || . 

In view of the unexpected nature of this compression effect, and its possible repercussions, as well as the extreme 
complexity of strong field general relativistic hydrodynamics, it is of course imperative that there be an independent 
confirmation of the existence of neutron star compression before one can be convinced of its operation in binary 
' systems. This is particularly important since many authors |]5|- |llj| have searched for but not observed this effect in 
Q\ , Newtonian tidal forces ||, first post-Newtonian (1PN) dynamics J6|,^0|, tidal expansions ]7]||, or in binaries in which 
rigid corotation has been imposed ||[llj . In ref. Q it has been argued that none of the above works could or should 
q-i have observed the effect, since the compression only dominates over tidal forces for stars with a realistic compaction 
ratio that are in a nearly irrotational hydrodynamic state (i.e. little spin relative to a distant observer). Indeed, 
irrotational stars may be a likely configuration near the final orbits as corotation would demand an unrealistically 
large viscosity in neutron stars 113] . 

With this in mind, it is of particular interest that a new formalism 13 Dj] has been proposed in which the hydrostatic 



quasi-equilibrium of irrotational stars in a binary can be solved independently of the complexities of (3+1) numerical 
' relativistic hydrodynamics. Thus, this is the first opportunity to independently test the hydrodynamic result. In this 
paper we report on results from an application of this independent formalism to binary neutron stars. We show that 
the interior density can increase as irrotational stars approach. 

Recently Flanagan |H[ ] has pointed out an inconsistency in the solution of the shift vector in previous calculations 
This problem is not present in the following simulations. More recently, a paper by Bonazzola, Gourgoulhon 
and Marck p9| presented a similar calculation, using stars with a smaller compaction ratio than the one in the present 
work. In their work the central density seemed to go unchanged through out the entire spiral motion until the very 
last orbits, in which it decreases slightly. Note, that though when this behavior does not significantly manifest the 
compression effect, neither does it obey any of the decreasing trends predicted by approximations mentioned above 
(i.e. Newtonian tidal forces, first post-Newtonian (1PN) dynamics, and tidal expansions). This clearly suggests 
the danger of approximating binary neutron star systems in their last stable orbits by weak field and/or corotating 
schemes. The binary systems studied in the present paper are composed by identical stars with compaction ratio 
higher than the one from Bonazzola et al, making them more relativistic objects. 



II. THE MODEL 



A full discussion of the CFC method can be found in [||. We use the (3+1) spacetime slicing as defined in the 
Arnowitt-Deser-Misner (ADM) formalism pf]|,pl|. Utilizing Cartesian x,y,z isotropic coordinates, proper distance is 
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expressed 



ds 2 = -(a 2 - (3 n f3 n )dt 2 + 2(3 n dx n dt + lns dx n dx 8 



(1) 



where the lapse function a describes the differential lapse of proper time between two hypersurfaces. The quantity j3i 
is the shift vector denoting the shift in space-like coordinates between hypersurfaces and 7^ is the spatial three-metric. 
The Latin indices go from 1 to 3. 

Using York's (3+1) formalism the initial value equations can be written as follows; the Hamiltonian constraint 
equation can be written 



R = 1671-p + K ns K ns - K 2 , 

where R is the Ricci scalar curvature, K ns is the extrinsic curvature, and p is the mass-energy density. 
The momentum constraint have the form p2| 

D n (K ni - j ni K) = 8wS i , 

where D n is the three-space covariant derivative and S l is derived from the stress-energy tensor 

where n v is a normal vector to a spatial slice. 

The CFC method restricts the spatial metric 7^ to the form 
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(2) 



(3) 



(4) 



(5) 



where the conformal factor </> is a positive scalar function. This approximation simplifies greatly the equations. It is 
motivated in part by the fact that the gravitational radiation in most systems studied so far is small compared to 
the total gravitational mass. The CFC is a frequently employed approach to the initial value problem in numerical 
relativity. Its application here is consistent with the quasi-equilibrium orbit approximation. 

The CFC leads to a set of elliptic equations for the metric components. Using Eq. (||) in combination with the 
maximal slicing condition tr(K) = 0, we get the following equations for (f> and (cap)'- 



V 2 </> = -4wpi 



(6) 



V 2 (a0) = 4ttp2 , 

where the Vi represent flat-space derivatives and the source terms are 



Pi 
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P0 W 2 + p e[T{W 2 - 1) + 1] + —K ns K 

1D7T 



(7) 



(8) 
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Pv(ZW 2 - 2) + p Q e[3T(W 2 + 1) - 5] 



16tt 



(9) 



where po is the rest-mass density, e the internal energy per unit of rest mass, T the adiabatic index, and W a 
generalization of the special relativistic 7 factor |Q| . A solution of Eq. (Q) determines the lapse function after Eq. (j|) 
is used to determine the conformal factor. 
The shift vector /3 1 can be decomposed fl23|: 



This is introduced into Eq. (||) to obtain the following two elliptic equations: 

V 2 x = V„B" , 



(10) 



(11) 
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V 2 B l = 2V n \n(a(j)- e )K m - 16na(/) 4 S i 



(12) 



An equation for the extrinsic curvature K % i is derived using the time evolution equation and the maximal slicing 
condition 

^' = ^(V^ + V^-^V n /3 n ) , (13) 

where K ij = 4> W K^. 

We assume that the matter behaves like a perfect fluid with a stress-energy tensor 

T^ = (p (l + e) + P) u^ + Pg^ , (14) 

and use a polytropic equation of state (EOS) 

P = kp T , (15) 

with P the pressure and k a constant taken as 1.13 x 10 5 erg cm 3 g -2 . This gives a maximum neutron-star gravitational 
mass of 1.46 Mq. In these simulations we consider two equal-mass neutron stars with a gravitational mass of 
mc = 1.42Mq each in isolation. This corresponds to a baryon mass of raj = 1.55 Mq. For the grid resolution of this 
study (~ 40 zones across the star in average) we obtain a central density in isolation of p c = 1.75 x 10 14 g cm" 3 . The 
compaction ratio for these stars is (mc/R)oo — 0.19, similar to the one of the stars considered in ||?]. As pointed 
out in Q| it is important to study realistically compact neutron stars. Otherwise Newtonian tidal forces can dominate 
over the relativistic effects one desires to probe. 



III. IRROTATIONAL STARS 



The method we use to determine the internal structure of stars in irrotational quasi-equilibrium configurations is 
essentially that originally proposed by Bonazzola, Gourgoulhon & Marck |l3| and as simplified by Teukolsky p| . 
The derivation of the relevant equations can be found in those papers. The essential ingredient of this approach is 
that, if the fluid vorticity is zero in the inertial frame, the specific momentum density per baryon can be written as 
the gradient of a scalar potential, 

hu^ = V M V > (16) 

where it M is the covariant four velocity and h is the relativistic enthalpy, h = 1 + e + j^. The potential ip can be 
obtained from the solution of a Poisson-like equation: 



D i D i i, = D i ^-- (D^-^B^Diln^j , (17) 

where Di are spatial covariant derivatives, n the baryon number density, and 

A = C + B , fl J ^a[/i 2 + (J}>Ai/')] 1/2 , (18) 

where C is a constant and B z is the shift vector in the rotating frame, B % — (3 l + (ui x r)% where u> is the angular 
velocity of orbital motion. 

Equation ([17]) must be solved by imposing a boundary condition at the stellar surface, 



a 2 ' 



= . (19) 

surf 



For irrotational stars the Bernoulli integral for the matter distribution then becomes: 



A 2 

h 2 = _(DV)AV'+-2 • (20) 



Equation ( pfj| ) uniquely determines the equilibrium structure of the stars. 

We also compute stars in constant corotation. In this case the relativistic Bernoulli equation can be written [p|JTl|l 
h/u° = constant. 

Solutions are obtained for specific values of the coordinate distance between stars and the total baryonic mass. This 
permitted us to construct a constant baryonic-mass sequence of orbits with a minimum number of code runs. This 
sequence is a collection of semi-stable orbits that are connected by the inspiral motion of the stars. 
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IV. NUMERICAL ALGORITHMS 



The problem consists essentially in the numerical solution of a set of elliptic equations. This set of equations is 
solved numerically using an iterative algorithm based upon an specially designed elliptic solver. This method consists 
of a combination of multigrid algorithms and domain decomposition techniques |ll| , p4[ ] and utilizes a code which was 
developed independently of the hydrodynamics code of [|l]j^] . 

The cases of corotating jll| and irrotational systems share the same set of equations ([|refalpha|ll], and |lj) for 
the metric fields. The irrotational systems, however, demand the solution of an extra elliptic equation for the 
description of the stellar internal structure. This poses a very special problem since the boundary condition ( |l9f ) is 
to be satisfied on the stellar surface and not on the grid boundaries as for the rest of the elliptic equations. This is 
particularly difficult to implement numerically since the stellar surface is a spheroid embedded in a Cartesian grid. 
We approach this problem in a way that resembles the traditional method for solving elliptic equations found in many 
textbooks (e.g. pq] ). Equation ( fL7| ) can be written as 

V 2 i/; = -47rp 3 , (21) 



where the source term p% is derived from the rhs of (0j]). The gradient of the baryon number density evaluated at the 



stellar surface in equation (|19|) is proportional to the unit vector s orthogonal to the same surface (i.e. Din oc — I). 

surf 

To simpli fy t he problem, we approximate the unit vector s by the radial unit vector r [ p6[ . Using this approximation, 
equation (|19|) reduces to 



dip 
Or 



= gm > (22) 

surf 



where g is a function of the position vector on the surface r. 

We next write ip in the form ip — ipp + ipL such that ipp satisfies the Poisson equation (|2l| ) with homogeneous 
boundary conditions on the outer limits of the numerical grid and ipL satisfies a Laplace equation with the following 
boundary conditions at the stellar surface 



dipi 



Or 



, (23) 

V Or J surf 



surf V dr 



Since ipL is a solution to a Laplace equation, it can be written as a linear combination of the spherical harmonics Yj m 

^L = J2Y.( Aimrl+Bimr ~ {l+1) ) Yi ™ (2 4 ) 

Z=0 m=-l ^ ' 

Only the terms up to I < 4 are included since numerical tests show that higher order multipoles are negligible. The 
coefficients Bi m must vanish to avoid singularities at the coordinate origin and the term A$o can be ignored since it 
only adds a constant to a potential. Thus, from equation (G4) we obtain 



dip 



^EE lAlmr 1 - 1 ^ . (25) 



1=1 m=-l 

If we consider a spherical surface S of coordinate radius r s , we can use the orthogonality condition of spherical 
harmonics to derive an equation for the Ai m coefficients 

J dUY lm ^ 



A lm = 7-|=i / dil Ytm-^r ■ (26) 



The numerical algorithm solves first the Poisson equation ( pi| ) for ipp using the same solver employed for the metric 
fields elliptic equations. Then dipp/dr is evaluated numerically on the spherical surface that best fits the actual stellar 
surface. This field is now used in combination with (|24|) and (^6|) to obtain ipL, and thus ip. 

This method has been tested on a Newtonian binary system where the stellar structure is described by a spherically 
symmetric mass density profile such that an analytical solution for the corresponding Poisson equation is known. A 
similar test was conducted with the density distribution obtained from one of the solutions for the irrotational cases 
presented in figure § The average relative difference between the corresponding benchmark solutions (the analytical 
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solution in the first case and the numerical solution obtained using the metric fields elliptic solver for the second case) 
and the solution provided by the ip solver was always less than 0.1%. 

Several other tests were performed to study the numerical stability of the ip algorithm and the dependence of the 
results on the boundary conditions ([l9|). In one test, the rhs of equation ( ^2|) was multiplied by scaling factors 0.8 and 
1.2. Another test included the re-scaling of ipL by these same factors. Such perturbations make only a small change 
in the results. The fractional change in central density (Ap/ p^^) deviates from the unperturbed result by less than 
0.02 for all cases. 

Our calculated quantities for the corotating sequence also agree well with values from Baumgarte et al. 
Differences in the fractional change in central density Ap/p 1 ^ 1 ^ are of the order of 10% [^7j. We attribute this 
difference to an effect of numerical resolution. 

Finally, the code was tested against the Newtonian irrotational sequences obtained by Uryu and Eriguchi |28| and 
Lai et al. p9[ . The code was stripped down of all the relativistic terms to reduce it to the Newtonian limit. The 
comparison was done for orbits with the same separation distance between stellar centers, d. The resulting values for 
the separation distance, the total energy, total angular momentum, and orbital angular frequency (and the relative 
difference with respect to HH|) were ^ = 3.53, & = -1-1419 (0.1%), J = 1.278 (4%), and fi = 0.2408 (3%). 



V. RESULTS 



We have found solutions to the initial value equations for semi-stable circular orbits for a binary system of identical 
neutron stars extending from the post-Newtonian regime to the innermost orbit for which we obtain a stable solution 
(probably near the last stable orbit fTl]] ). Figure [l] shows the change Ap in the central rest-mass density relative to 
the central density Pq of an isolated star. Results are plotted as a function of the proper distance between stellar 
centers. The numerical error for the irrotational points was conservatively estimated to be < ±0.01 based upon the 
convergence tests performed on the code, while for the corotating ones it was estimated to be < ±10% based upon 
the comparison with the results from For both sets of points numerical results are consistent (within numerical 
error) with the expectation ||f7]j|] that the change in central density should be small and negative in the Newtonian 
limit. 
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FIG. 1. Change in central density relative to a single isolated star (Ap/pg"^) as a function of the proper distance between 
stellar centers. 
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Figure [l] shows the main result of this paper. There is a clear unambiguous difference between irrotational and 
corotating stars. In both cases the central density approaches the isolated-star limit at large orbital separations. For 
closer orbits and higher frequency, the central density decreases for corotating stars. This is consistent with the results 
of Refs. |Q,0,O. The exact opposite is true, however, for irrotational stars. We have checked that the same trend 
emerges in a plot of the average density, so this effect could not be an artifact of the stellar center being a special 
point. 

Although the qualitative result of increasing density as the orbit shrinks is reproduced here, we note that the 
magnitude of the effect is less than that observed in the previous unconstrained hydrodynamic calculations ^,|). 
Also, these results do not reproduce the (l/£) 6 dependence of central density with separation as expected in a tidal 
expansion analysis This is no surprise, since this analysis is for stars of relatively large ratio of neutron star 

radius to orbital separation (R/L). In the present work, (R/L) goes up to 0.5, thus falling outside the range of validity 
of a tidal expansion. 

In a recent paper Bonazzola et al. |19| performed a similar calculation using different numerical methods. They 



report no changes in the central density for irrotational binaries for most of the sequence and a slight decrease at the 
end, when the stars approach their closest separation distance. This calculation was done for a binary with stars with 
compaction ratio of 0.14, which represent stars more extended than the ones portrayed here. Preliminary simulations 
performed by our group for a similar binary system show results compatible with those of Bonazzola et al. and will 
be reported in an upcoming paper |jo"| . 

The calculations reported here also exhibit a trend of decreasing binding energy with increased orbital frequency 
(and decreasing separation) as expected. The orbital frequencies we derive for both irrotational and corotating stars 
are very close to the Newtonian frequency. Figure || shows the binding energy as a function of the proper distance. 
Note that the sequence does not present a minimum which characterizes the point of dynamical instability of the 
system. 
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FIG. 2. Binding energy of the system, denned as the half the gravitational mass of the binary minus the gravitational mass 
of one star in isolation, divided by the latter mass (Mo — Maa ) /M^ as a function of the proper distance between stellar centers 
in km. 



In summary, the present independent study has confirmed the qualitative result of increasing central density as an 
irrotational binary orbit decays, for stars with compaction ratio of 0.19. What is still required, however, is a fully 
dynamical relativistic hydrodynamics simulation. This awaits the completion of the neutron-star Grand Challenge 
project. 



G 



VI. ACKNOWLEDGMENTS 



Work at University of Notre Dame supported by NSF grant PHY-97-22086. Work at the Lawrence Livermore 
National Laboratory performed in part under the auspices of the U. S. Department of Energy under contract W-7405- 
ENG-48 and NSF grant PHY-9401636. 



VII. REFERENCES 



[1] J.R. Wilson and G.J. Mathews, Phys. Rev. Lett. 75, 4161 (1995). 

[2] J.R. Wilson, G.J. Mathews, and P. Marronetti, Phys. Rev. D54, 1317 (1996). 

[3] G. J. Mathews and J. R. Wilson, Astrophys. J. 482, 929 (1997). 

[4] G J. Mathews, P. Marronetti and J. R. Wilson, Phys. Rev. D58, 043003 (1998). 

[5] D. Lai, Phys. Rev. Lett. 76, 4878 (1996). 

[6] A. G. Wiseman, Phys. Rev. Lett. 79, 1189 (1997); M. Shibata, Prog. Theor. Phys. 96, 317 (1996); M. Shibata, K. 
Taniguchi, and T. Nakamura, Prog. Theor. Phys. 128, 295 (1997); K. Taniguchi and T. Nakamura, Prog. Theor. Phys. 96, 
693 (1996); D. Lai and A. G Wiseman, Phys. Rev. D54, 3958 (1996); W. Ogawaguchi and Y. Kojima, Prog. Theor. Phys. 
96, 901 (1996). J. C. Lombardi, F. A. Rasio, and S. L. Shapiro, Phys. Rev. D56, 3416 (1997); P. Brady and S. Hughes, 
Phys. Rev. Lett. 79, 1186 (1997). 
E. Flanagan, Phys. Rev. D58 (1998) 124030. 
K. Thorne, Phys. Rev. D58 (1998) 124031, 

T. W. Baumgarte, G B. Cook, M. A. Scheel, S. L. Shapiro, and S. A. Teukolsky, Phys. Rev. Lett. 79, 1182 (1997); Phys. 
Rev. D 57, 6181 (1998), Phys. Rev. D 57, 7299 (1998). 

M. Shibata, T. W. Baumgarte, and S. L. Shapiro, Phys. Rev. D58 (1998) 023002. 
P. Marronetti, G J. Mathews, and J. R. Wilson, Phys. Rev. D58 (1998) 042822. 
L. Bildsten and C. Cutler Astrophys. J. 400, 175 (1992). 

S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Phys. Rev. D5 6, 7740 (1997). 
S. A. Teukolsky, Astrophys. J. 504, 442 (1998), ( |gr-qc/9803082| ). 
M. Shibata, Phys. Rev. D58 (1998) 024012. 
H. Asada, Phys. Rev. D57 (1998) 7292. 



E. Gourgoulhon, Phys. Rev. D, submitted (1998), (^r-qc/9804054j) . 
E. Flanagan, Phys. Rev. Lett. 82 1354 (1999). 

S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Phys. Rev. Lett. 82 892 (1999). 

R. Arnowitt, S. Deser, and C. W. Misner, in Gravitation, edited by L. Witten (Wiley, New York, 1962) p. 227. 

J. W. York, Jr., in Sources of Gravitational Radiation, edited by L . Smarr (Cambridge Univ. Press, Cambridge, England, 

1979), p. 83. 

C. R. Evans, Ph.D. Thesis, University of Texas, 1985; P. Anninos, D. Hobill, E. Seidel, and L. Smarr, Phys. Rev. Lett. 71, 
2851 (1993). 

J. M. Bowen and J. W. York Jr., Phys. Rev. D21, 2047 (1980). 

P. Marronetti and G. J. Mathews, J. Comp. Phys, submitted (1998) . 

J. D. Jackson, Classical Electrodynamics, (John Wiley, New York, 1963). 

This is consistent with the fact that the solutions obtained for irrotational stars show very little distortion from a spherical 
shape. 

This comparison was done on the limited range of overlap between our calculations and those of Ref. || . 
K. Uryu and Y. Eriguchi, Mon. Not. Roy. Astron. Soc. 296 (1998) LI. 

D. Lai, F. A. Rasio, and S. L. Shapiro, Astrophys. J. 420, 811 (1994). 
P. Marronetti, G. J. Mathews, and J. R. Wilson. In preparation. 



7 



